Generated from rmarkdown file: “code/behav/_behav.rmd“. This rmd file executes ‘child’ R-scripts within this directory. These child scripts contain analysis code. Their output is captured and ‘knitted’ back into this report. To run the code, you can either execute this .rmd file (i.e., ‘knit’ to html), or source each child script individually. (These child scripts should be configured to run on their own when sourced directly.)

1 analysis-set subjs: model construction, estimation of stroop effects

1.1 RT analysis

1.1.1 plot

plot(behav$rt)
abline(h = 250)
abline(h = 3000)  ## marks beginning of subsequent trial

behav$is.implausible.rt <- behav$rt > 3000 | behav$rt < 250

grid.arrange(
  
  behav %>%
    
    filter(acc == 1, !is.implausible.rt) %>%
    
    ggplot(aes(rt)) +
    geom_density(fill = "slateblue", alpha = 0.3) +
    
    labs(title = "all correct trials") +
    theme(legend.position = "none"),
  
  behav %>%
    
    filter(acc == 1, !is.implausible.rt) %>%
    
    ggplot(aes(rt)) +
    geom_density(aes(fill = trial.type), alpha = 0.3) +
    
    labs(title = "by congruency") +
    scale_fill_brewer() +
    scale_color_brewer() +
    theme(legend.position = c(0.5, 0.5)),

  ncol = 2
  
)

behav %>%
  
  filter(acc == 1, !is.implausible.rt) %>%
  
  full_join(group_by(., subj) %>% summarize(r2norm = qqr2(rt)), by = "subj") %>%
  
  ggplot(aes(sample = rt)) +
  stat_qq(alpha = 0.8, size = 1) +
  stat_qq_line(size = 0.25) +
  
  facet_wrap(vars(subj)) +
  geom_text(
    aes(
      x = -1.25, y = 2000,
      label = paste0("r^2 = ", round(r2norm, 3))
    ), size = 3, color = "grey50"
  ) +
  theme_minimal(base_size = 10) +
  theme(axis.title = element_blank()) +
  labs(
    title    = "QQ rt versus normal",
    subtitle = paste0("overall r^2 with normal = ", round(qqr2(behav$rt), 3))
  )

849971 and 161832 have odd patterns in RT distribution. Strong clustering at 500 ms. This is highly consistent with an artifact we found in other RT data from the same microphone/preprocessing method.

1.1.2 handling likely-artifactual RTs

behav %>%
  
  filter(acc == 1, !is.implausible.rt) %>%
  
  ggplot(aes(x = rt, group = subj)) +
  geom_density(data = . %>% filter(!subj %in% c("849971", "161832")), color = rgb(0, 0, 0, 0.15), size = 2) +
  geom_density(data = . %>% filter(subj %in% c("849971", "161832")), color = "firebrick1", size = 2) +
  labs(
    title    = "subjects with bimodal distributions (likely artifactual) in red"
  )

behav$is.artifactual.rt <- FALSE
behav$is.artifactual.rt[behav$subj %in% c("849971", "161832") & behav$rt < 500] <- TRUE

behav %>%

  filter(acc == 1, !is.implausible.rt, subj %in% c("849971", "161832")) %>%
  full_join(group_by(., subj) %>% summarize(r2norm = qqr2(rt)), by = "subj") %>%
  
  ggplot(aes(sample = rt)) +
  stat_qq(alpha = 0.8, size = 0.4) +
  stat_qq_line(size = 0.25) +
  geom_hline(yintercept = 500) +
  facet_wrap(vars(subj)) +
  theme_minimal(base_size = 10)

1.1.3 fit RT model

## make sure to exclude validation set!
behav.rt.aset <- behav %>% filter(acc == 1, !is.implausible.rt, !is.artifactual.rt, is.analysis.group)

fit1.het <- lme(
  rt ~ trial.type, 
  random  = ~ trial.type | subj,
  data    = behav.rt.aset,
  weights = varIdent(form = ~ 1 | subj),
  control = cl1,
  method  = "REML"
)

## define "far" outliers as points with resids more extreme than 3 IQR

behav.rt.aset$resid.p <- resid(fit1.het, type = "p")
behav.rt.aset$is.far.out <- farout(behav.rt.aset$resid.p)

## exclude far outliers and refit model

fit1.het.trim <- update(fit1.het, data = behav.rt.aset %>% filter(!is.far.out))
fit1.het.trim.ml <- update(fit1.het.trim, method = "ML")  ## for model comparisons

1.1.4 test heterogeneity of level-I variance

fit1.hom.trim.ml <- lme(
  rt ~ trial.type, 
  random  = ~ trial.type | subj,
  data    = behav.rt.aset %>% filter(!is.far.out),
  control = cl1,
  method  = "ML"
)

(rt.hom.v.het <- anova(fit1.hom.trim.ml, fit1.het.trim.ml))
##                  Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## fit1.hom.trim.ml     1  6 132710.8 132754.1 -66349.38                        
## fit1.het.trim.ml     2 54 128508.2 128898.3 -64200.10 1 vs 2 4298.551  <.0001

1.1.5 test variance in stroop effect

fit0.het.trim.ml <- update(fit1.het.trim.ml, random  = ~ 1 | subj)
(rt.stroopvar <- anova(fit0.het.trim.ml, fit1.het.trim.ml))
##                  Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## fit0.het.trim.ml     1 52 128591.7 128967.4 -64243.87                        
## fit1.het.trim.ml     2 54 128508.2 128898.3 -64200.10 1 vs 2 87.53103  <.0001

1.1.6 estimate cross-run reliability

fit1.het.run.trim <- lme(
  rt ~ trial.type * run, 
  random  = ~ trial.type * run | subj,
  data    = behav.rt.aset %>% filter(!is.far.out),
  weights = varIdent(form = ~ 1 | subj),
  control = cl1,
  method  = "REML"
)
summary(fit1.het.run.trim)
## Linear mixed-effects model fit by REML
##  Data: behav.rt.aset %>% filter(!is.far.out) 
##        AIC      BIC    logLik
##   128353.5 128808.6 -64113.76
## 
## Random effects:
##  Formula: ~trial.type * run | subj
##  Structure: General positive-definite, Log-Cholesky parametrization
##                     StdDev    Corr                
## (Intercept)         117.38521 (Intr) trl.ty run   
## trial.typeincon      33.46107  0.008              
## run                  42.03440  0.315 -0.007       
## trial.typeincon:run  11.09470 -0.430  0.084  0.708
## Residual             82.02205                     
## 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | subj 
##  Parameter estimates:
##      107321      123117      130114      132017      138837      141422 
##   1.0000000   2.2042468   1.1374682   1.9329308   1.8295447   1.7368152 
##      150423      158136      160830      161832      165032      171330 
##   2.0390034   1.0513835   4.5850017   5.1846430   2.6913126   1.8198504 
##      178243      178950      182436      197449      203418      204319 
##   2.8558033   1.3055077   3.2531529   2.1067865   1.1163342   2.5114102 
##      300618      317332      346945      352738      448347      580650 
##   1.0329272   1.1376803   1.4761791   2.5163718   1.6665064   3.3350494 
##      594156      601127      765864      814649      843151      849971 
##   1.0921172   1.9365139   1.0036256   1.7891142   2.1310886   2.7519015 
##      873968      877168 DMCC1971064 DMCC2442951 DMCC3062542 DMCC4854075 
##   1.2534443   2.1359043   1.1221483   0.7311733   1.7021922   1.2050789 
## DMCC5009144 DMCC6418065 DMCC6484785 DMCC6661074 DMCC6671683 DMCC6721369 
##   0.8460941   0.9366789   1.7051477   1.3783664   0.9051912   0.9766868 
## DMCC7297690 DMCC7921988 DMCC8033964 DMCC8050964 DMCC9441378 DMCC9478705 
##   1.3092629   2.0403024   0.8997665   1.6099094   2.4952706   1.5766959 
## DMCC9953810 
##   1.5094231 
## Fixed effects: rt ~ trial.type * run 
##                        Value Std.Error    DF  t-value p-value
## (Intercept)         787.5277 18.064731 10089 43.59476  0.0000
## trial.typeincon      75.4953  9.052337 10089  8.33986  0.0000
## run                  22.3310  7.293891 10089  3.06160  0.0022
## trial.typeincon:run   2.9946  5.016553 10089  0.59695  0.5506
##  Correlation: 
##                     (Intr) trl.ty run   
## trial.typeincon     -0.226              
## run                  0.047  0.318       
## trial.typeincon:run  0.103 -0.725 -0.200
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -3.8309061 -0.6703863 -0.1488609  0.4980128  4.6330611 
## 
## Number of Observations: 10141
## Number of Groups: 49
## get unconditional / marginal covariances

m <- rbind(c(0, 1, 0, 0), c(0, 1, 0, 1))  ## contrast matrix
tau.trim <- getVarCov(fit1.het.run.trim)
(tau.trim <- m %*% tau.trim %*% t(m))  ## covariance matrix
##          [,1]     [,2]
## [1,] 1119.643 1150.952
## [2,] 1150.952 1305.353
(r.marginal.trim <- cov2cor(tau.trim)[1, 2])  ## correlation
## [1] 0.952036
svd(getVarCov(fit1.het.trim), nv = 0L)$d  ## covmat not degenerate (eigvals > 0)
## [1] 22654.493  1429.647

1.2 error analysis

1.2.1 test variance in stroop effect

behav <- behav %>% mutate(error = 1 - acc)

fit.error0 <- glmer(
  error ~ trial.type + (1 | subj),
  behav %>% filter(response.final != "unintelligible", is.analysis.group), 
  family  = binomial,
  control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1E9))
)
fit.error1 <- glmer(
  error ~ trial.type + (trial.type | subj), 
  behav %>% filter(response.final != "unintelligible", is.analysis.group), 
  family = binomial,
  control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1E9))
)
(er.stroopvar <- anova(fit.error0, fit.error1))
## Data: behav %>% filter(response.final != "unintelligible", is.analysis.group)
## Models:
## fit.error0: error ~ trial.type + (1 | subj)
## fit.error1: error ~ trial.type + (trial.type | subj)
##            Df    AIC    BIC  logLik deviance  Chisq Chi Df Pr(>Chisq)
## fit.error0  3 1729.8 1751.6 -861.89   1723.8                         
## fit.error1  5 1733.0 1769.3 -861.52   1723.0 0.7416      2     0.6902

1.2.2 estimate cross-run reliability

fit.error1.run <- glmer(
  error ~ trial.type * run + (trial.type * run | subj),
  behav %>% filter(response.final != "unintelligible", is.analysis.group),
  family = binomial,
  control = glmerControl(optimizer = "Nelder_Mead", optCtrl = list(maxfun = 1E9))
)
## Warning in checkConv(attr(opt, "derivs"), opt$par, ctrl = control$checkConv, :
## Model failed to converge with max|grad| = 0.047928 (tol = 0.001, component 1)
summary(fit.error1.run)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: error ~ trial.type * run + (trial.type * run | subj)
##    Data: 
## behav %>% filter(response.final != "unintelligible", is.analysis.group)
## Control: 
## glmerControl(optimizer = "Nelder_Mead", optCtrl = list(maxfun = 1e+09))
## 
##      AIC      BIC   logLik deviance df.resid 
##   1740.3   1841.9   -856.1   1712.3    10516 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.4186 -0.1392 -0.0843 -0.0528 19.4999 
## 
## Random effects:
##  Groups Name                Variance Std.Dev. Corr             
##  subj   (Intercept)         0.3950   0.6285                    
##         trial.typeincon     1.9923   1.4115    0.75            
##         run                 0.2797   0.5288    0.85  0.43      
##         trial.typeincon:run 0.8166   0.9036   -0.81 -0.95 -0.63
## Number of obs: 10530, groups:  subj, 49
## 
## Fixed effects:
##                     Estimate Std. Error z value Pr(>|z|)    
## (Intercept)          -5.2767     1.0806  -4.883 1.04e-06 ***
## trial.typeincon      -0.1373     1.2254  -0.112    0.911    
## run                  -0.2723     0.7203  -0.378    0.705    
## trial.typeincon:run   0.8742     0.7929   1.102    0.270    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) trl.ty run   
## tril.typncn -0.858              
## run         -0.923  0.823       
## trl.typncn:  0.828 -0.943 -0.904
## convergence code: 0
## Model failed to converge with max|grad| = 0.047928 (tol = 0.001, component 1)
plot(rePCA(fit.error1.run)$subj)

1.3 extract blups and draw figures

## get RT blups
blups <- as.data.frame(coef(fit1.het.trim))
blups %<>% rename(congr = "(Intercept)", stroop = "trial.typeincon") %>% tibble::rownames_to_column("subj")

## bind with error blups

blups %<>%
  
  full_join(
    
    data.frame(
      subj = rownames(coef(fit.error1)$subj),
      er.logit.stroop = coef(fit.error1)$subj$trial.typei,  ## extract logits
      er.logit.congr  = coef(fit.error1)$subj[["(Intercept)"]]
    ) %>%
      mutate(
        er.logit.incon = er.logit.stroop + er.logit.congr,  ## logit of error on incon trials
        ##  blup stroop effect in units percent error:
        stroop.er = (logit2prob(er.logit.incon) - logit2prob(er.logit.congr)) * 100
      ) %>%
      dplyr::select(subj, stroop.er),
    
    by = "subj"
    
  )
## Warning: Column `subj` joining character vector and factor, coercing into
## character vector
## draw figure

plot.behav <- 
  
  blups %>%
  
  mutate(subj = factor(subj, levels = subj[order(stroop, decreasing = TRUE)])) %>%
  dplyr::select(subj, stroop, stroop.er) %>%
  reshape2::melt() %>%
  
  
  filter(!is.na(value)) %>%
  
  ggplot(aes(subj, value)) +
  facet_grid(
    rows = vars(variable), scales = "free", switch = "y",
    labeller = as_labeller(c(stroop = "Response time", stroop.er = "% error"))
  ) +
  
  geom_segment(aes(x = subj, y = 0, xend = subj, yend = value), color = "grey50", size = geom.line.size) +
  geom_point(fill = "grey30", color = "black", shape = 21, size = geom.point.size) +
  coord_capped_cart(left = "both") +
  
  xlab("Subject") +
  theme(
    panel.grid       = element_blank(), 
    panel.border     = element_blank(),
    panel.background = element_blank(),
    strip.placement  = "outside",
    strip.background = element_blank(),
    strip.text       = element_text(size = axis.title.size),
    axis.line.y     = element_line(size = axis.line.size*0.6),
    axis.text.y     = element_text(size = axis.text.size),
    axis.text.x     = element_blank(),
    axis.ticks.y    = element_line(size = axis.line.size),
    axis.ticks.x    = element_blank(),
    axis.title      = element_text(size = axis.title.size),
    axis.title.y = element_blank()
  )
## Using subj as id variables
ggsave(here("out", "behav", "stroop-blups.pdf"), height = 2.5, width = figwidth)

## write

write.csv(blups, here("out", "behav", "stroop_blups_rt_group201902.csv"),  row.names = FALSE)
write.csv(behav, here("out", "behav", "behavior-and-events_group201902_with-subset-cols.csv"),  row.names = FALSE)

behav.mod.objs <- list(
  
  fit1.het.trim = fit1.het.trim,
  er.stroopvar = er.stroopvar,
  rt.stroopvar = rt.stroopvar,
  rt.hom.v.het = rt.hom.v.het,
  r.marginal.trim = r.marginal.trim
  
)
saveRDS(behav.mod.objs, here("out", "behav", "mod_objs.RDS"))

2 microphone comparison

2.1 read data

2.2 preliminary look

behav %>%
  
  mutate(subj = factor(as.factor(subj), levels = micinfo[order(micinfo$mic), "subj"])) %>%
  
  ggplot(aes(subj, rt, color = mic)) +
  geom_point(position = position_jitter(width = 0.1), size = 0.5, alpha = 0.3) +
  geom_boxplot(position = position_nudge(1/3), notch = TRUE, width = 0.2) +
  
  scale_color_brewer(type = "qual", palette = 2) +
  
  theme(
    axis.text.x = element_text(angle = 90, hjust = 0, vjust = 0),
    legend.position = c(0.8, 0.8)
  )
## Warning: Removed 110 rows containing non-finite values (stat_boxplot).
## Warning: Removed 110 rows containing missing values (geom_point).

behav %>%
  
  group_by(mic, subj) %>%
  summarize(freq = sum(rt == 0, na.rm = TRUE)) %>%
  
  ggplot(aes(mic, freq)) +
  geom_point(size = 2) +
  
  labs(y = "frequency of rt == 0 per subj")

behav %>%
  
  group_by(mic, subj, error) %>%
  summarize(freq = sum(error, na.rm = TRUE)) %>%
  
  ggplot(aes(mic, freq)) +
  geom_point(size = 2, position = position_jitter(width = 0.1), alpha = 0.4) +
  
  labs(y = "frequency of errors per subj")

behav %>%
  group_by(mic, subj, acc.final) %>%
  summarize(freq = n()) %>%
  ggplot(aes(mic, freq)) +
  facet_grid(cols = vars(acc.final)) +
  geom_point(size = 2, position = position_jitter(width = 0.1), alpha = 0.4) +
  labs(y = "frequency of acc.final codings per subj")

behav %>%
  
  ggplot(aes(sample = rt, color = mic)) +
  
  stat_qq(alpha = 0.3, size = 0.15) +
  stat_qq_line(size = 0.25) +
  
  facet_wrap(vars(subj)) +
  scale_color_brewer(type = "qual", palette = 2)
## Warning: Removed 110 rows containing non-finite values (stat_qq).
## Warning: Removed 110 rows containing non-finite values (stat_qq_line).

behav.mean.sd <- behav %>%
  
  group_by(mic, subj, session) %>%
  
  summarize(rt.mean = mean(rt, na.rm = TRUE), rt.sd = sd(rt, na.rm = TRUE))

grid.arrange(
  
  behav.mean.sd %>%
    ggplot(aes(mic, rt.mean)) +
    facet_grid(vars(session)) +
    geom_point(position = position_jitter(width = 0.1)) +
    geom_boxplot(position = position_nudge(1/3), notch = FALSE, width = 0.2),
  
  behav.mean.sd %>%
    ggplot(aes(mic, rt.sd)) +
    facet_grid(vars(session)) +
    geom_point(position = position_jitter(width = 0.1)) +
    geom_boxplot(position = position_nudge(1/3), notch = FALSE, width = 0.2)
  
)

2.3 model

2.3.1 differences in mean stroop effect btw mics?

  • leave in validation-set data for additional power
behav.rt.mic <- behav %>% filter(!is.implausible.rt, !is.artifactual.rt, acc == 1, mic != "unknown")

fit <- lmer(
  rt ~ trial.type * mic + (trial.type | subj),
  behav.rt.mic
)
summary(fit)
## Linear mixed model fit by REML. t-tests use Satterthwaite's method [
## lmerModLmerTest]
## Formula: rt ~ trial.type * mic + (trial.type | subj)
##    Data: behav.rt.mic
## 
## REML criterion at convergence: 156411.5
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -3.1439 -0.5426 -0.1402  0.3388 10.7989 
## 
## Random effects:
##  Groups   Name            Variance Std.Dev. Corr
##  subj     (Intercept)     12992    113.98       
##           trial.typeincon  1233     35.11   0.29
##  Residual                 27684    166.38       
## Number of obs: 11948, groups:  subj, 57
## 
## Fixed effects:
##                          Estimate Std. Error     df t value Pr(>|t|)    
## (Intercept)                718.81      33.39  54.93  21.528  < 2e-16 ***
## trial.typeincon             91.14      12.29  53.67   7.413 9.09e-10 ***
## micmicro                   104.34      37.58  54.95   2.776   0.0075 ** 
## trial.typeincon:micmicro   -12.70      13.85  53.84  -0.917   0.3633    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr) trl.ty micmcr
## tril.typncn  0.153              
## micmicro    -0.888 -0.136       
## trl.typncn: -0.136 -0.888  0.153
fit.het <- lme(
  rt ~ trial.type * mic, 
  random  = ~ trial.type | subj,
  weights = varIdent(form = ~ 1 | subj),
  data    = behav.rt.mic,
  control = cl1
)
summary(fit.het)
## Linear mixed-effects model fit by REML
##  Data: behav.rt.mic 
##        AIC      BIC    logLik
##   152783.6 153256.4 -76327.81
## 
## Random effects:
##  Formula: ~trial.type | subj
##  Structure: General positive-definite, Log-Cholesky parametrization
##                 StdDev    Corr  
## (Intercept)     114.40479 (Intr)
## trial.typeincon  32.88231 0.289 
## Residual        158.98579       
## 
## Variance function:
##  Structure: Different standard deviations per stratum
##  Formula: ~1 | subj 
##  Parameter estimates:
##      115825      123117      130114      130518      132017      138837 
##   1.0000000   1.3680243   0.9491809   0.7273882   1.1413267   0.9937305 
##      141422      150423      158136      165032      171330      178243 
##   1.0374857   1.3434624   0.6237094   1.6795750   0.9487168   1.5708969 
##      178647      178950      197449      203418      204319      233326 
##   1.1122717   0.7281854   1.2069523   0.7358722   1.6734130   2.0684924 
##      300618      317332      346945      352738      393550      448347 
##   0.5321997   0.6442652   0.8708645   1.6949101   0.7676516   1.1962199 
##      580650      594156      601127      765864      814649      843151 
##   2.0987369   0.6359272   1.1860198   0.5632435   0.9543656   1.2030346 
##      849971      873968      877168 DMCC1328342 DMCC1596165 DMCC1971064 
##   1.5060506   0.7877884   1.3172820   0.3994237   1.6274925   0.8102444 
## DMCC2442951 DMCC2803654 DMCC2834766 DMCC3062542 DMCC4191255 DMCC4854075 
##   0.4411251   0.6736451   0.6987490   0.9810488   0.6522941   0.6798883 
## DMCC5009144 DMCC5195268 DMCC5775387 DMCC6371570 DMCC6418065 DMCC6661074 
##   0.5090925   0.5882508   1.4021924   0.9881116   0.5390242   0.7096957 
## DMCC6671683 DMCC6705371 DMCC6721369 DMCC7297690 DMCC8033964 DMCC8050964 
##   0.5231536   0.7451889   0.5573742   0.7791344   0.4987172   0.8518369 
## DMCC9441378 DMCC9478705 DMCC9953810 
##   1.3638986   0.8554003   0.9115365 
## Fixed effects: rt ~ trial.type * mic 
##                             Value Std.Error    DF   t-value p-value
## (Intercept)              721.5696  33.38987 11889 21.610437  0.0000
## trial.typeincon           86.4733  11.16612 11889  7.744258  0.0000
## micmicro                 100.3302  37.59825    55  2.668481  0.0100
## trial.typeincon:micmicro  -7.0790  12.63690 11889 -0.560188  0.5754
##  Correlation: 
##                          (Intr) trl.ty micmcr
## trial.typeincon           0.182              
## micmicro                 -0.888 -0.162       
## trial.typeincon:micmicro -0.161 -0.884  0.178
## 
## Standardized Within-Group Residuals:
##        Min         Q1        Med         Q3        Max 
## -4.8864061 -0.6227512 -0.1716170  0.4302583 10.0354679 
## 
## Number of Observations: 11948
## Number of Groups: 57
  • fomri microphone recorded faster RTs
  • no observed impact of microphone on size of mean stroop effect ### differences in variance between mics?
fit.het.ml  <- update(fit.het, . ~ ., method = "ML")
fit.het.mic.ml <- update(fit.het.ml, . ~ ., weights = varIdent(form = ~ 1 | mic))
fit.hom.mic.ml <- update(fit.het.ml, . ~ ., weights = NULL)
(anova.mic <- anova(fit.het.ml, fit.het.mic.ml, fit.hom.mic.ml))
##                Model df      AIC      BIC    logLik   Test  L.Ratio p-value
## fit.het.ml         1 64 152811.9 153284.8 -76341.96                        
## fit.het.mic.ml     2  9 156422.1 156488.6 -78202.06 1 vs 2 3720.196  <.0001
## fit.hom.mic.ml     3  8 156456.2 156515.3 -78220.10 2 vs 3   36.081  <.0001
modobjs <- list(
  anova.mic = anova.mic,
  mic.model.means = fit.het
)

saveRDS(modobjs, here("out", "behav", "mod_objs_mic.RDS"))
  • microoptics microphone recorded more variable RTs

3 validation-set subjs: estimation of stroop effects

3.1 read data

3.2 model, plot, and write

## initial fit

behav.vset.rt <- behav %>% filter(acc == 1, !is.na(rt), rt < 3000, rt > 250, !is.analysis.group)

fit.vset <- lme(
  rt ~ trial.type, 
  random  = ~ trial.type | subj,
  data    = behav.vset.rt,
  weights = varIdent(form = ~ 1 | subj),
  control = lmeControl(maxIter = 1e5, msMaxIter = 1e5, niterEM = 1e5, msMaxEval = 1e5),
  method  = "REML"
)

## trim and re-fit

behav.vset.rt$resid.p <- resid(fit.vset, type = "p")
behav.vset.rt$is.far.out <- farout(behav.vset.rt$resid.p)
fit.vset.trim <- update(fit.vset, subset = !is.far.out)

## model error

fit.error.vset <- glmer(
  error ~ trial.type + (trial.type | subj), 
  behav %>% filter(response.final != "unintelligible", !is.analysis.group) %>% mutate(error = 1 - acc), 
  family = binomial,
  control = glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1E9))
)
## boundary (singular) fit: see ?isSingular
summary(fit.error.vset)
## Generalized linear mixed model fit by maximum likelihood (Laplace
##   Approximation) [glmerMod]
##  Family: binomial  ( logit )
## Formula: error ~ trial.type + (trial.type | subj)
##    Data: 
## behav %>% filter(response.final != "unintelligible", !is.analysis.group) %>%  
##     mutate(error = 1 - acc)
## Control: glmerControl(optimizer = "bobyqa", optCtrl = list(maxfun = 1e+09))
## 
##      AIC      BIC   logLik deviance df.resid 
##    680.3    711.4   -335.2    670.3     3661 
## 
## Scaled residuals: 
##     Min      1Q  Median      3Q     Max 
## -0.4709 -0.1620 -0.1031 -0.0076  9.6958 
## 
## Random effects:
##  Groups Name            Variance Std.Dev. Corr 
##  subj   (Intercept)     17.87    4.227         
##         trial.typeincon 10.12    3.181    -1.00
## Number of obs: 3666, groups:  subj, 17
## 
## Fixed effects:
##                 Estimate Std. Error z value Pr(>|z|)  
## (Intercept)      -11.961      4.715  -2.537   0.0112 *
## trial.typeincon    7.983      4.657   1.714   0.0865 .
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Correlation of Fixed Effects:
##             (Intr)
## tril.typncn -0.998
## convergence code: 0
## boundary (singular) fit: see ?isSingular
## extract predictions

blups.vset <- as.data.frame(coef(fit.vset.trim))
blups.vset %<>% rename(congr = "(Intercept)", stroop = "trial.typeincon") %>% tibble::rownames_to_column("subj")

## bind with error blups

blups.vset %<>%
  
  full_join(
    
    data.frame(
      subj = rownames(coef(fit.error.vset)$subj),
      er.logit.stroop = coef(fit.error.vset)$subj$trial.typei,  ## extract logits
      er.logit.congr  = coef(fit.error.vset)$subj[["(Intercept)"]]
    ) %>%
      mutate(
        er.logit.incon = er.logit.stroop + er.logit.congr,  ## logit of error on incon trials
        ##  blup stroop effect in units percent error:
        stroop.er = (logit2prob(er.logit.incon) - logit2prob(er.logit.congr)) * 100
      ) %>%
      dplyr::select(subj, stroop.er),
    
    by = "subj"
    
  )
## Warning: Column `subj` joining character vector and factor, coercing into
## character vector
## plot

plot.behav.vset <- 
  
  blups.vset %>%
  
  mutate(subj = factor(subj, levels = subj[order(stroop, decreasing = TRUE)])) %>%
  dplyr::select(subj, stroop, stroop.er) %>%
  reshape2::melt() %>%
  
  filter(!is.na(value)) %>%
  
  ggplot(aes(subj, value)) +
  facet_grid(
    rows = vars(variable), scales = "free", switch = "y",
    labeller = as_labeller(c(stroop = "Response time", stroop.er = "% error"))
  ) +
  
  geom_segment(aes(x = subj, y = 0, xend = subj, yend = value), color = "grey50", size = geom.line.size) +
  geom_point(fill = "grey30", color = "black", shape = 21, size = geom.point.size) +
  coord_capped_cart(left = "both") +
  
  xlab("Subject") +
  theme(
    panel.grid       = element_blank(), 
    panel.border     = element_blank(),
    panel.background = element_blank(),
    strip.placement  = "outside",
    strip.background = element_blank(),
    strip.text       = element_text(size = axis.title.size),
    axis.line.y     = element_line(size = axis.line.size*0.6),
    axis.text.y     = element_text(size = axis.text.size),
    axis.text.x     = element_blank(),
    axis.ticks.y    = element_line(size = axis.line.size),
    axis.ticks.x    = element_blank(),
    axis.title      = element_text(size = axis.title.size),
    axis.title.y = element_blank()
  )
## Using subj as id variables
plot.behav.vset

## write and save

ggsave(here("out", "behav", "stroop-blups-validation.pdf"), height = 2.5, width = figwidth/2)
write.csv(blups.vset, here("out", "behav", "stroop_blups_rt_group201902_validation.csv"),  row.names = FALSE)
saveRDS(fit.vset.trim, here("out", "behav", "fit1-het-trim_group201902_validation_set.RDS"))